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different functional forms for the internal spin and overall phase of the order parameter. 
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^ = 1. For ^ > 2 we obtain topological as well as non-topological solutions defined 
by the asymptotic radial dependence. For arbitrary values of £ the non-topological 
solutions include bright ring-vortices which explicitly demonstrate the confining effects 
of the Dirac operator. We arrive at solutions through an asymptotic Bessel series, 
algebraic closed-forms, and using standard numerical shooting methods. By including a 
harmonic potential to simulate a finite trap we compute the discrete spectra associated 
with radially quantized modes. We demonstrate the continuous spectral mapping 
between the vortex and free particle limits for all of our solutions. 
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Vortices appear in physical settings which span a wide range of energy scales and 
disciplines ranging from hurricane phenomena [1] and turbulent flow in classical 
fluids [2], to superfluidity in ^He [Hj and general quantum fluids lais], down to the 
smallest length scales of Bogomol’nyi-Prasad-Sommerfleld states in supersymmetric 
field theories [HI. Vortices are relevant from a technological standpoint in quantum 
computing for example [7j, as well as in more theoretical areas of research such as 
galactic halos in Bose-Einstein condensate (BEC) theories of dark matter jH]. In BECs 
vortices as stable rotating solutions of the nonlinear Schrodinger equation (NLSE) are 
ubiquitous P [To] . Indeed, the presence of a persistent quantized rotation may be 
considered a defining property of superfluidity Although a variety of types of vortices 
are possible in atomic condenstates m, the Laplacian inherent to the Schrodinger 
Hamiltonian constrains the number of observable vortex structures in the BEC. For 
instance, a BEC made up of spin-F bosons gives rise to a {2F + l)-component order 
parameter which nevertheless solves a multi-component NLSE. An alternative method 
of constructing a multicomponent BEC is to add the internal degrees of freedom by 
placing a condensate in a two-dimensional honeycomb optical lattice [lain]. Atoms 
are condensed into the lowest energy Bloch state then translated to a corner of the 
Brillouin zone using laser assisted Bragg scattering [T3| . At wavelengths large compared 
to the lattice constant, the microscopic details of the lattice manifest as an additional 
SL(2, C) spin group symmetry, i.e., a Dirac point structure emerges [T2] - 

In a previous paper we derived a BEC version of the nonlinear Dirac equation 
(NLDE) for the case of weak interparticle interactions [I2| . This has attracted attention 
from diverse fields of research [lailSliniilHlilglEolEllEaEalElEslEniEZlEHlI^ 
EOlEIlISlISSlESISSlEniEZlEHlISnilinillll. Solutions of the NLDE are effectively long 
wavelength limit lattice envelopes, and so can cover a large number of sites. The healing 
length in our effectively quasi-2D system is typically about 10 times the lattice constant, 
and we require the healing length to be small compared to the size of the cylindrical 
container. In typical experiments, the number of lattice sites is on the order of 100 in a 
linear direction, thus vortex solutions of the NLDE may be described accurately. In this 
article we study relativistic quantum vortices in the superffuid phase of a Bose gas at the 
Dirac point of a honeycomb optical lattice focusing on explicit vortex solutions of the 
NLDE. Our vortices are solutions of equations reminiscent of relativistic systems, but 
pertain here to ultracold atomic gases, where we work in the mean-field limit throughout. 
Solution types differ by asymptotic conditions on the amplitude for radius r much greater 
than the healing length, far from a vortex core, or r much less than the healing length 
deep inside a vortex core. A general schematic for vortices in our problem is shown in 
Fig.Il Solutions may also differ by use of the internal degree of freedom in the spinor 
structure and with respect to quantization of rotation. However, we find that spinor 
components within a particular solution always differ by one unit of rotation. 

Quantized rotation in a (2 -|- l)-dimensional spin-1/2 order parameter is best 
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Figure 1. (color online) Relativistic vortex in the honeycomb lattice. The vortex 
current is depicted by red arrows with the core towards the interior of the red circle. 
The relativistic current is quantized around the core and incorporates the overall U(l) 
phase into the internal SU(2) rotation. 


characterized in terms of the topology of the circulation or phase winding around a 
singular core. Consider an overall macroscopic phase which wraps around the circular 
boundary consisting of radial lengths much larger than the healing length. When the 
internal spin degrees of freedom remain topologically unconstrained we simply referred 
to this as a vortex. In contrast, when the wrapping around the circular boundary 
includes the internal spin degrees of freedom in a nontrivial way the order parameter is 
usually called a texture or skyrmion min]. Throughout our work we adhere to the 
case where the Berry phase is spliced into the overall U(l) symmetry upon condensation 
of the order parameter. Thus, we do not consider solutions with a macroscopic yr-phase 
disinclination. We categorized our solutions according to their topological properties 
based on symmetries of the order parameter manifold. These symmetries must take 
into account the internal spin degrees of freedom as well as the overall U(l) phase. 
In particular, localized solutions may asymptotically approach either a local minima, 
maxima, or saddle point of the effective potential in the mean-field energy functional. 
The first leads to topological categorization, whereas the latter two cases are not 
topologically protected configurations. 

The presence of the Dirac operator in the NLDE adds new features over conventional 
spinor BECs. In some cases the derivative terms act to defocus the solution leading to 
conventional “dark” vortices with density minima in an otherwise uniform background. 
In other cases the cross-mixing of spinor components in the Dirac term has a confining 
effect resulting in a “bright” vortex made of spinning localized density rings over a zero 
density background. These localized solutions are similar to ones found in ordinary 
BECs with attractive interactions [131 mi, although we emphasize that in the present 
work we consider only repulsive contact interactions. It is significant that our solutions 
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Figure 2. (color online) Schematic of relativistic nonlinear Dirac vortices. At the 
Dirac point (horizontal direction) changing the interaction, lattice constant, particle 
hopping, and density tunes the order parameter between bright vortices (red dots 
over yellow) with different quantized rotation i G Z+. The peaks and valleys 
denote topological transitions between elements of the fundamental group 7ri(5'^) 
around the vortex core. Changing the chemical potential away from the Dirac point 
(vertical direction) transitions the order parameter into skyrmion solutions with a 
single quantum of rotation i = 1. Higher energy then leads to multiple phase winding 
states i >2 (red dots over blue), the energy increasing with winding i. 


apply to certain spin-orbit coupled BECs [45] . In the low-energy limit of these theories 
the linear (Dirac) kinetic term is dominant and one may safely neglect the second 
order (Schrodinger) contribution. Tuning the interactions to eliminate cross terms in 
the spinor components leads to our NLDE. The various types of NLDE vortices are 
summarized in Fig[^ 

To put relativistic vortices in a more general context, the large symmetry 
groups encountered in high energy physics and cosmology admit vortex solutions as 
fundamental ingredients. Examples include cosmic snperstrings [IHl HU HE], mirror 
symmetry of Calabi-Yan manifolds m [50], and brane world scenarios EH Vortices 
are generic to gauge theories where spontaneous symmetry breaking occurs. This 
includes both the Abelian case @81152], as with Nielsen-Olesen vortices in the Abelian 
Higgs model, in addition to the non-Abelian case [S3], non-Abelian Yang-Mills 
for example. Supersymmetric held theories which exhibit weak-strong duality rely 
fundamentally on the presence of solitons and vortices, as these provide the necessary 
degrees of freedom to make the dualities possible jHll HH] . 

Our results are organized as follows. In Sec. we express the full time-dependent 
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NLDE in (2+l)D in radial form appropriate for stationary vortices. This reduces 
the problem to the solution of two coupled first order nonlinear ordinary differential 
equations. In Sec. symmetries of the order parameter are determined in order to 
classify the various types of solutions. Section develops the Lagrangian formulation 
of the nonlinear Dirac equation and connects to the standard theory of spontaneous 
symmetry breaking. In Sec. we use the mean-field energy functional to examine 
in detail the interplay between the effective potential and Dirac kinetic terms to gain 
deeper insight into quasi relativistic vortices. In Sec. we solve the NLDE analytically 
using asymptotic Bessel and algebraic methods. Density and phase plots depicting each 
solution type are provided. Section connects our solutions to vortices in the nonlinear 
Schrodinger equation with correction terms, through a semiclassical reduction. In Sec.|^ 
we present numerical methods and solutions for vortices. This approach allows us to 
explore the full vortex landscape by tuning the various physical parameters using the 
method of numerical shooting. Section [^presents vortex solutions and spectra in a weak 
harmonic trapping potential building on results from the previous sections. Finally, in 
Sec. [ 10 ] we conclude. 


2. Vortices in the nonlinear Dirac equation 


In this section we express the NLDE in the appropriate coordinates for finding vortex 
solutions and analyze the solution space at the energy functional level. We show that 
in different regions of solution space the effective potentials switch between focusing 
and defocusing forms, providing conditions for either bright (localized ring) or dark 
(conventional) vortices. The NLDE describes the dynamics of a four-spinor in (2 -|- 1)- 
dimensions, = (d^+, d'-)^, with the upper (-I-) and lower (—) two-spinors relating to 
opposite K and K' points of the honeycomb lattice and ^ C^. The full NLDE 

can be expressed as 

[ ihYd^, -UY^ 1 ^ = 0 , ( 1 ) 

\ ^= 0,1 / 

where the 4x4 interaction matrices are given by 

M, = t(7»)^+i(-irW^7V, (2) 

constructed to give the correct cubic nonlinearites, local to each spinor component m. 
The matrices {p, = 0, 1, 2) are the usual Dirac matrices in the chiral 

representation [55]. Equation Q describes a gapless theory which corresponds to 
massless interacting Dirac spinors. Since the interactions do not couple different spinor 
components, Eq. Q can be split into two sets of equations (one for each of the K and 
K' points) of the form 


ihdt^+ = [ -ihci a-V + UY j T+ , 


(3) 


i=0,l 
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where we use the Pauli matrices a = {a^-, (Jy, cr^), and is the nonzero 2x2 

submatrix along the diagonal of M, given explicitly by 

Mf« = . (4) 

Note the presence of the effective speed of light q and interaction strength U. These 
qnantities are defined in terms of the more fundamental lattice and atomic parameters 
as Q = tha-\/3/2h and U = LzQfV 3-s/3a^/8, where a and t^ are the lattice spacing and 
hopping energy, and L^, g, h are the vertical oscillator length, two-body interaction, and 
average three-dimensional atomic density, respectively. The associated long-wavelength 
healing length associated with the NLDE is then found to be ^Dirac = tha\/Z/2U. For a 
complete list of the relevant renormalized physical parameters with experimental values 
see [Hj. 

We seek planar stationary solutions to Eq. (|^, and thus we require the upper 
two-spinor to be expressed in factorized form d'+(r,t) = exp{—ifj,t/h) ['0^(r), 'Vls(r)] , 
where g is the chemical potential of the system. The spinor components t^A(B)(r) here 
are complex functions on the plane, i.e., C. For stationary solutions in 

rectangular coordinates, Eq. (§ further reduces to 

- ihci (d^ - idy) y) + U I'lpAix, y)f ipAix-, y) = gi^Aix-, y) , (5) 

-ihci (9^ -F idy) i^Aix, y) + U I'lpsix, y)f ipsix., y) = gipsix-, y ), (6) 

with the full solution expressed as a linear combination of solutions from Dirac points 
K and Kk To obtain cylindrically symmetric solutions with arbitrary integer phase 
winding, i.e. vortices, we transform to plane-polar coordinates. The time-independent 
spinor wavefunctions are then written as 

Mr ,«) = ± i , Mr, 0) = e^feir ), (7) 

where £ is the angular momentum quantum number (integral phase winding) and the 
radial functions fA(B) are real functions of r. The asymmetric dependence on the winding 
nnmber i between spinor components in Eq. Q is reqnired. In polar coordinates, the 
derivatives in Eq. (§-(@ acquire factors of exp(—id) and exp(fd), respectively. An offset 
factor of exp(—z6*) in 'ij^A is then reqnired in order to separate the rotational and radial 
dependence. Applying the method of separation of variables forces the particular form 
of angular dependence in Eq. 0 . Equations. 0-0 then reduce to the two coupled 
radial equations 

- hci + D |/n(r)|Vn(r) = p/n(r), (8) 

hci (dr + ^^^fA(r) + U \fBir)f fsir) = gfB{r) • (9) 

Equations (|^-@ comprise the key system of conpled nonlinear equations of this article. 
In the sections that follow we will solve these using combinations of analytical and 
nnmerical techniques for various asymptotic and core bonndary conditions. Localized 
solutions of Eqs. @-([^ may be categorized by their asymptotic forms, i.e., whether the 
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amplitudes /a(b) have zero or nonzero limits far from the core, a relevant distinction 
in our problem. In particular, we note that in conventional BECs with repulsive 
interactions the nonlinearity is always defocnsing which leads to a vortex with the 
familiar constant ambient profile that vanishes only at the core. This contrasts radically 
to the physics in our problem, which we will see allows for (bipartite) vortices with 
nonzero asymptotic profiles in addition to ones with asymptotically vanishing profiles. 
The latter structures are the bright vortices or ring-vortices previously mentioned. 


3. Symmetries of the order parameter 


To gain a foundational understanding of vortex solutions of the NLDE, we study first 
the symmetries of the order parameter. Condensed bosons near the Dirac point of a 
honeycomb lattice are characterized by U(l) symmetry breaking in the full many body 
theory, and transfer from the gronnd state to the metastable state at the edge of the 
Brillouin zone edge mi. The order parameter is comprised of 2 complex amplitudes 
'ifj = '4’b)'^ associated with the bipartite lattice, forming a representation of the 

two-dimensional special unitary gronp SU(2) = U(1)0 ® Spin(2). The Spin gronp in 
2D forms a double cover of the rotation group SO(2) = U(l). Spin components xpA 
and tl^B are conpled throngh the kinetic part of the single particle Hamiltonian bnt are 
independent with respect to onsite interactions. Thus, the full symmetry group G for a 
Dirac-point BEC that leaves the mean-field energy fnnctional invariant is 


G'^U(1)^0U(1)s0(Z2)5 , 


( 10 ) 


where the subscripts 0 and S denote the gauge and spin degrees of freedom, respectively. 
The 2D representation of G is generated by the 2x2 unit and Pauli matrices 
1, a^, (Ty, whose Lie group is given by exp {—i(l>) x exp {—in ■ d/2), where a = {a^, (Ty), 
h = (cosd, sind), a G R, d is the polar angle, and (f) is the gange degree of freedom. In 
matrix form, G can be parameterized as 


G ^ 


0 

gie/2 



( 11 ) 


where again 0 is the gauge angle, 9 is the polar angle, and where the discrete cyclic 
(off-diagonal matrix) part of G is isolated by setting 0 = 0 = 0. One sees that an 
arbitrary order parameter 0, obtained by applying G to a representative spinor 00 , has 
an additional discrete symmetry H. This is evident by taking 0 —?• 0-|-7r and 9 ^ 6 + 27r 
in Eq. (11) described by the discrete symmetry (Z 2 )^^, leading us to deduce the order 


parame 


)er manifold M = G/H 


For the gapless system considered here, we generally reqnire both sublattices to 
be occupied, so that a representative standard spinor is given by 0o = (1, 1)^- An 
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arbitrary order parameter tl; is obtained by applying the symmetry transformation G to 
■ 00 ) which gives 


_ j{<t>-ei2) 


ijj = e 



= e 


iV 


Ae 


(13) 


Requiring the order parameter to be single valued under spatial rotations places a 
constraint on the gauge angle. In particular, we require (p to contain a part that depends 
on the polar angle by defining p' = p-\- 0/2 and taking p' as the new gauge parameter. 


Note that we could have factored Eq. ((T3|) another way and obtained 


-,-10 


1 


= e*'^' 


-,-i9 


(14) 


in which case we would define p' = p -\- 9/2. The discrete choice of redefinition 
p' = p 9/2 splits the order parameter space into right and left chiral spinor 


manifolds displayed in Eq. (13) and Eq. (14), respectively. Imposing the condition of 


single valuedeness for the order parameter has effectively eliminated the extra discrete 
symmetry Einally, computing the fundamental group of M then yields the 

direct sum 7ri(M) = 7ri(G') = 0 l- 2 ,s- 


4. Lagrangian formulation and spontaneous symmetry breaking 

So far we have focused primarily on the symmetries of the order parameter associated 
with the Dirac point structure. We now examine the effect of the nonlinearity associated 
with contact atomic interactions. Before we obtain exact vortex solutions to the 
nonlinear Dirac equation it is helpful to study the theory at the Lagrangian level in 
order to gain insight into the interplay of kinetic and interaction terms and to acquire a 
qualitative sense of the kinds of vortex solutions one should expect. From a field theory 
viewpoint, Eq. 0 comprises the Euler-Lagrange equations to the Lagrangian density 

£ = {p^dtpA + P^dtpB) - ihn2v>ci [p*A{d^ - idy)'pB + 'P^ii.dx + idy)'pA] 

(15) 

where Eq. 0 is obtained through the standard prescription 
dC \ dL 


d„ 


dp 


+ 




dp* 


= 0 . 


(16) 


Note that a factor of the two-dimensional average density n 2 D appears in Eq. (15), which 


gives C the correct units of energy (see Ref. [SS])- The positive sign on the second term 
is correct, since we are taking derivatives with respect to the conjugate of the field and 
its derivatives. The more common approach is to differentiate the Lagrange density with 
respect to the field and its derivatives and take the conjugate of the resulting equations 


afterwards. This leads to the same result. Eq.(15) describes the dynamics of two self¬ 


interacting, scalar fields coupled through the spatial derivative terms, with interaction 
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% = - C , 


(17) 


where it a and ttb are the canonical momenta associated with xpA and ips respectively, 
defined in the usual way by 
dC 


= 


’ 


(18) 


The Lagrangian density in plane polar coordinates reads 

1 


C = ihn2D {'lp*Adt'lpA + - ihn2DCi 


U 


il}*Ae"\dr 


d0)'lpB + + i-d0)'lpA 




(19) 


Since we consider only the case of repulsive interaction where U > 0, and given the 
linear derivatives in the Lagrangian density, the condensate can lower its energy by 
finding favorable combinations of fields and gradients that render the derivative terms 
negative. This leads to the formation of a relativistic vortex on an otherwise uniform 
positive energy background. Spontaneous symmetry breaking is exhibited in Eq. (19) 
by expressing the field energy as the sum of a uniform contribution equal to the chemical 
potential of the system p, plus an additional piece that accounts for the presence of a 
vortex. Substituting the ansatz t) = t) and t) = e t) 

in Eq. (19), where the factor of 2 accounts for an even sublattice asymptotic particle 
distribution, we obtain 

£ = ihn2G{v*AdtVA + v^dtUs) 


ihn2GCi 


* 


V aC 


{dr 


i - d0)vB + v*Be + i - de)vA 




U. 


(Kr+Kn-y(Kr+Kr). 


Reading off the effective potential energy density in Eq. 


we find 


Te = V\vA? + \VB\^)-^(\vA\‘ + \VB\'). 


U . 


( 20 ) 


( 21 ) 


Decomposing the sublattice order parameters in terms of amplitude and phase va(b) = 
^Pa(b) we see that the minima of the effective potential energy are located at 

four points in spin space with infinite U(l) degeneracy at each point: 



Thus, we see that topological solutions are constrained beyond the symmetry 
classifications presented so far. At ultra low temperatures the nonlinearity arising 
from local onsite atomic interactions locks the order parameter into a one of four 
possible spin configurations characterized by different relative internal phases. However, 
taking onto account the phases 4>a{b)) the symmetry group of the NLDE nonlinearity 
reduces to U(l),^^ 011 ( 1 ) 0 ^ 0 ~^a,Bi where the discrete transformation simply exchanges 
va and vb in Eq. (|^. So far, our analysis shows that topological solutions of the 










Vortex solutions and spectra in a weak harmonic trap 


10 


NLDE are possible and should involve mapping of one or both U(1 )^(b) phases to the 
circle at spatial infinity. When both phases are involved, the vortex core describes 
tunneling between potential minima (n™™, n™") ^ —(n™", In contrast, when 


only a single U(l) mapping is involved, a vortex core describes a single tunneling 
(n™, n™") (—n™, n™") or (n™", —n™). We note that non topological solutions are 

also possible. These are solutions that extremize Eq. (15), but asymptotically approach 


saddle points or local maxima of Eq. (21). Such solutions are possible even in the 


presence of strictly repulsive interactions. This is because of the rich structure contained 


in the spin-orbit type (kinetic) coupling in Eq. (15), the effects of which we shall examine 
in further detail. 


5. Energy functional analysis for vortices 


To study the interplay between the effective potential and kinetic terms of the NLDE, 
we would like to map out the energy landscape for a vortex configuration in more detail. 
Here, we work with the mean-field energy functional for a generic vortex by eliminating 
the angular dependence and incorporating centripetal terms into the effective potential 
energy. Assuming the vortex form in Eq. 0 and taking ^ = 0 in our analysis, the total 
energy is given by 


where, 

Kfr[/A,/n] 


E = f dr [hn 2 r>ci {-IbIa + fAfs) + Ws] , 




^2U 


1 ) + 


IbIa 




(23) 


(24) 


and we regulate the energy by introducing an upper cutoff R. Note that radial 
derivatives are indicated by prime notation. Evidently the fAf's fsfA terms 
provide an attractive force for configurations where sgn(/^/g) > 0 and sgn(/^/B) < 0, 
respectively. In the /^, /^-plane, the centripetal term in V^fi dominates for small r and 
forms a saddle point at the origin (of the /^, /^-plane) becoming singular when r = 0. 
Because of the saddle-point, points on the and fs axes have zero potential energy 
but are unstable. If we want to study solutions that begin at the saddle-point, i.e., 
/n(0) = /b(0) = 0, or on the /^-axis, say, near the saddle-point, we must include the 
full contribution from the kinetic terms. In fact, we see that a path along the /^-axis 
defined by p(t):= tes, 0 < t < tf, will change the total energy by AE = Qt//^(0). 
This can be made arbitrarily large and negative by adjusting the value of at r = 0. 
The I4ff-landscape flattens rapidly as r increases from zero, so that it may be possible 
to attempt a solution for which /a(0) = 0, /n(0) > 0, f^iO) < 0. 

If we consider finite energy solutions at large r, and must go to zero and 
the centripetal term may be neglected. In this limit, we find that Veff has nine local 
extrema, four of which are minima for the total energy E and thus associated with 
topological solutions. We obtain the following asymptotic (large r) critical points for 
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in the /^j/s-plane: , \/2pL/U) - local minimum; {y^2fi/U, 0) - saddle 

point; {^J2|l/U, —^J2|l/U) - local minimum; {—^/2|J/U, ^J2jl/U) - local minimum; 
{—^2ii/U, 0) - saddle point; —^2ii/U) - local minimum; (0, y'2n/U) - 

saddle point; (0, —y^%I/U) - saddle point; (0, 0) - local maximum. For our chosen 
case study with £ = 0, i.e., with vorticity in only of spinor component, we then expect 
topological solutions with properties /a(0) = 0, |/^(0)| ^ 0, |/n(0)| > ^/2pLJU^ fsi^) — 
0, and lim^-j-ood/ul, |/s|) = (\/2/i/{7, -y/2/i/{7). This behavior describes a Mermin-Ho 
vortex. Note that more general topological solutions exist for arbitrary winding with 
rotation in both spinor components which must both vanish at r = 0, in order to 
maintain the finite energy requirement. Moreover, we find that several non-topological 
localized solutions exist as well, with asymptotic limits coinciding with local extrema 
that do not minimize the potential energy. Such solutions may satisfy the conditions 
/u(0) = 0, |/^(0)| 0, /b( 0) 0, /^(O) = 0, and lim^^oolU,/ b) = (0,0), which 

describes a ring-vortex/soliton. Another possible non-topological solutions occurs for the 
same initial conditions but with asymptotic behavior lim^^-ood/Al,/ b) = (\/2/r/17, 0), 
which describes the vortex/soliton solution. 

The appearance of non-topological solutions that are nevertheless localized is 
a consequence of the various possible combinations of configurations for the spinor 
components at the core of a vortex. Some of these result in attractive contribution 
to the total energy leading to self-trapping of the order parameter, which is particularly 
dramatic in the case of ring-vortex solutions. A look at the energy functional for Eqs. (|^- 
© will help us gaiu iusight iuto this effect. The NLDE euergy fuuctioual associated 
with a vortex with arbitrary circulatiou i is 


E = dr 


hn2T>Ci ( f^fB - IbIa + IbIa 


2i-l 


+ -r 

+ q/A 




B 


(25) 


Taking a different viewpoint, we can define the component-specific effective potentials 
as 


U^S = fAfB + fBfA 


+ “T (/I + /b) = /a/b + U£ + Ua-a , 


r 4 
2 £ — 1 U 

Vff = “ IbIa + /b/a -- 1- -J (/b + Ia) = ~fBfA + Ui + Ua-a , 


(26) 

^ \J a '•> A/ — 1 “I 1 “ u—u ; ( 27 ) 

where we have introduced condensed notation for the angular momentum and atom- 
atom interactions Ua-a-, respectively. The energy functional Eq. (25) reduces to the 
form 


E = 


dr [fiSjDC, (J'Jb - IbIa) + Ui, + Ug] 


(28) 


The effective potentials encapsulate the angular momentum, binary interactions, 

aud derivative-auiplitude biliuear terras. These effective poteutials exhibit trausitious 
betweeu focusiug aud defocusiug forms for differeut regious of spiuor solutiou space. 
A couveuieut way to characterize such trausitious is through the eigeuvalue of the 2D 
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Conditions for /a, fs and /y 


TJA 

JjB 

Dg 

1 - / a , /s, /y Pb > 0; IaPb >Ug + Ua-a] 

def. 

foe. 

null 

Ja > /y fB < Pb 

2. fA,fB,PB>0-,PA<0; /bIPaI 

) /a/s > Ug + Ua-a', 

foe. 

foe. 

neg. 

Ja > l/y; fB > Pb 

3. /a, /b, Pb > 0; /^ < 0; /bIPaI 

\, fAPB>Ug + Ua-a; 

foe. 

foe. 

neg. 

Ib ~ /a; \Pa\ ~ Pb 

4. /b? /a? /a ^ Oj /b ^ 0 


def. 

def. 

pos. 

5. fA, Pa> 0; fB, Pb< 0-, \fB\PA<Ui + Ua-a; 

def. 

def. 

pos. 


Table 1. Character of the effective potentials and the Dirac spectrum for different 
regions of the NLDE solution space. Effective potentials exhibit transitions 

between focusing (foe.) and defocusing (def.) forms depending on the values for 
spinor components and their derivatives. The character of D denotes the dominant 
contributing elements of the spectrum for given conditions on /a(b) and i.e., 

the largest contribution to a vortex solution may come from the positive (pos.) or 
negative (neg.) part of the spectrum or may be an equal admixture of the two (null). 
Note that conditions 1-5 do not exhaust all possibilities but are listed to demonstrate 
the variability needed to observe bright and dark vortices. 


massless radial Dirac operator defined as 


D, 


+ 0 j- 


(29) 


Denoting the continnous spectrum of Di as Di = Spec{Dg) C (—oo,+oo), we 
can further divide the spectrum into the negative, null, and positive subspaces 
{A_ I A_ G ( — oo , 0)}, {Ao I Ao = 0}, {A+ | A+ G ( 0 , +oo )}, respectively. The nature of 
the effective potentials and (and thus vortex solutions) is linked to the character 
of the subspace of Dg that dominates a particular region of the NLDE solution space. 
Table [T] displays the various regions based on analysis of the effective potentials in 


The conditions in Table [T] are not exhaustive but provide some evidence for 
the simultaneous presence of dark and bright vortices in the same physical system, 
particularly in the presence of repulsive atomic interactions. Upon obtaining explicit 
solutions, we will see that condition 1 and 2 describe the upswing and downswing 
behavior near the apex of the bright vortex and excitations of the asymptotically fiat 
vortex. This reversal in trend for the radial spinor profiles characterizes the self-trappin 
effect of focusing potentials. Condition 3 applies to all texture solutions whose spinor 
degrees of freedom are topologically constrained. Condition 4 gives rise to conventional 
dark vortices with arbitrary phase windings and non-vanishing tails but lacking the 
topological constraint of the skyrmion. 
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In this section we derive explicit vortex solutions of the NLDE using a combination 
of analytical and numerical techniques. Some of the solutions have been presented 
in our previous work [H] but without the mathematical detail in the present article. 
Different methods are applicable depending on whether the chemical potential is zero or 
finite-valued and if the vortex has a single unit of winding or arbitrarily large winding 
For arbitrary winding and nonzero chemical potential we use both an approximate 
asymptotic method and numerical shooting. In the case where £ = 0 or 1 exact analytical 
solutions are possible, and in particular when fi = 0 concise algebraic forms occur. 
For all of the solutions in this section we consider only the case of a spatially infinite 
condensate. In experiments, this is a good approximation when the healing length, and 
thus the core of the vortex, is small compared to the oscillator length associated with a 
trapping potential which sets the size of the BFC. 

6.1. Asymptotic Bessel solutions for large i 

Returning to the cylindrically symmetric form of the NLDF, Fqs. (§-(@) , we can see 
that for winding values £ > 2 both and /s must vanish at r = 0 due to the presence 
of the centrifugal terms. Thus, we will treat the special cases £ = 0,1 in a separate 
section. To obtain vortex solutions we require spatial derivatives to vanish at infinity; 
this means that for r —?• oo, Fqs. yield the asymptotic behavior 

lim fA{B){r) [|/n(s)(r)P - la/U] = 0 , (30) 

r—)-oo 

which implies 

lim \fA{B){r)\ =0, ±^Jii/U. (31) 

The important point to note here is whether the fA{B) tends to zero or a non-zero 
constant as r ^ oo. To determine asymptotic forms for the radial functions we first 
look for exponential decay to zero, or growth or decay to the constants i-yZ/r/t/, in 
which case we find no consistent asymptotic solution to Fqs. Algebraic decay to 

zero may be discerned by substituting /a(?’) — and /n(r) = r~^ into Fqs. @-([^ for 
a, (3 > 0, decay or growth to ±y//r/t/ by /a(?’) = r~'^±yl[i/U and /^(r) = r“^±y//r/t/. 
For the case of non-vanishing boundary conditions and £ > 2 the asymmetry in the 
angular momentum terms complicates the task of finding solutions. We recall that the 
factors dr + i/r and dr + {1 — i)/r act as index raising and lowering operators for the 
Bessel functions J„. One can verify that Bessel functions are nearly exact solutions in 
the case of a weak nonlinearity given by {7|d^p/fin2DQ 1. In this limit solutions of 
Fqs. (Il)-® are approximately Bessel functions of the first kind, Ji and J/_i. A weak 
nonlinearity only slightly modifies the Bessel form by perturbatively scattering J/ and 
Ji-i into higher £-valued states. 

However, the case of strong nonlinearity is different. In this case we expect the 
solution to deviate drastically from the Bessel form, especially for the non-vanishing 
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asymptotic profile of a vortex which interests us here. This heuristic discussion motivates 
a modified Bessel-Fourier expansion as our ansatz, which must include Bessel functions 
to all orders 


fA{f') = AF{r) e^^^A 


Q-o 


T ^ ^ ^nJn (^) 


n=l 


(32) 

Mr) = (B/A)Ur), (33) 

where A and B are normalization constants, Jn{r) is the Bessel function of the first 
kind of order n, ao and are the expansion coefficients. Anticipating an asymptotic 
expansion, we include the real function Q(r) in the complex prefactor in order to cancel 
oscillations in the Bessel series at large argument. In addition, we have included the 
radial function F(r), which we require to cancel the r~^ asymptotic decay of the Bessel 
functions. The series that we have chosen to use runs over the Bessel index rather than 
the usual form where the summation runs over the zeros of a single Bessel function 
with a fixed index. This choice of expansion is valid but does not offer the convenience 
of using the standard orthonormal relations for Bessel functions when computing the 
coefficients a„. 

Substituting the ansatz Eqs. (32)-(33) into Eq. (§-(0) , we then consolidate the 
angular momentum and derivative terms in the resulting series expansion by using 
the recurrence relations for Bessel functions: J^jr = (Jn-i + Jn+\)l^'^ and J!^ = 
{Jn-i — d„+i)/2. We obtain the recursion relations for the coefficients in addition to 
two first-order differential equation for the functions Q and F 

2 


F r hci 


ao + E OriJn 




hci oq -|- a^Jn 


iQ' = - 
jiA 


^ no(l - 
> + r 


^±Uaf? 

hci 


Oo 


T ^ ^ O'nJ'i 


Ti 


hciB 
fiB 


(1 - C)an = 


O-n+l f d F n + 1\ Un-l 


(1 - C)a„ 


2 


n -|- 1 
2-i + n 


+ 


H Uo T Ca^Jn 
hCl Oq F y E OjnJn 
£ — n F V 


F 


1 


n — 1 

2-£-n 


(34) 

(35) 

(36) 

(37) 


hciA 2 \ nFl J 2 \ n — 1 

In Eqs. (|M1)-(|^ we have absorbed a fraction 1 — C* < 1 from the chemical potential 
terms into the recursion relations. Solving the recursion relations leads to 

C fj, \ n ! 


A = F iB , 


'^V( 2 £-l) 


n 


n—1 


(38) 


Equations (34)-(35) are consistent only in the regime £ ^ 1. Moreover, since Q is a real 
function we require F'/F F ao£/r = 0. Choosing F(r) = r leads to Oq = —1/£. This 
choice of F effectively forces the solution to vanish at the origin and also cancels the 
r~^ behavior of the Bessel functions at long distances. For large £ and r ^ ^oirac in fhe 
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strong nonlinear regime pi/U = 1, Eqs. (|34|)-([l^ read 

|A|V2 


nci 




sinr 


^2n~ 


(^2n+l- 


cosr 


- 1 


hci 


i^r 

J^(-l)” (aansinr - a 2 n+iCosr) 

2 

- 1 


n 



which upon integration yields 

Q{r) 


, U 

hci 


(39) 

(40) 

(41) 


where after integration the oscillating part of Eq. ( |40| ) may be neglected compared to 

(42) 


the linear part. Note, that to obtain Eq. (39) we used Rayleigh’s formula 

'1 d ^ 


Ur) = (-r)” I 


smr 


Einally, the parameter C in Eq. (38) must be tuned to a critical value C'yoj-texi dehned 
so that the Bessel series converges and is nonzero for all values of i. We obtain 
Cyortex ~ (5/2) x ( 2 ^ — 1 ) with closed form solution given by 

n—1 


/^(r) = A{Ur/hci) 


n=l 


n ! 


n' 


— Jn{±Ur/hci). 


(43) 


The overall constant A is determined by normalizing the wavefunction to the particle 
number. We emphasize here that Eq. (43) is an approximate solution which applies only 
very near to or far from the vortex core, and for large winding number. Density and 
phase plots for the complex Bessel solution are shown in Fig. on the left with radial 
profile plots shown in Fig. |^a). 


6.2. Algebraic solutions 

Here we treat solutions for £ > 1 at zero chemical potential pi — 0. For zero chemical 
potential, algebraic closed forms exist for arbitrary values of the winding number. In 
contrast, for general values of p, closed forms exist only for I = 0 ,1, which we discuss 
in the last paragraph of this section. By setting /i = 0 in Eqs. @-0 we obtained 
exact non-topological vortex solutions. From a technical standpoint eliminating the 
chemical potential terms in the NLDE simplifies the problem considerably. This is 
because solutions to the homogeneous non-interacting equations are algebraic forms, 
not Bessel functions. Thus, vortex solutions do not connect to Bessel functions in 
the zero-interaction limit. We note here that the precise relationship between the 
chemical potential p and interaction strength U is determined by computing the 
spectrum of each solution confined within a harmonic trap. This is the subject of 
Sec. of the present article. To obtain algebraic solutions we start from an ansatz 
/a = + 'jr^y^'^ and likewise for fs with A ^ B, where the parameters to be 

determined are a, /?, 7 , 5 G R and A, B E C hj substitution into the NLDE. Substituting 
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Figure 3. (color online) Phase and density plots for the £ = 4 asymptotic Bessel 
topological vortex (left) and numerical topological vortex (right). (a,b) A-sublattice 
density and phase. (c,d) B-sublattice density and phase. 


this form for and fs into Eqs. ([^-(|^ and solving for A, B, a, (3 ,7 and S by matching 
coefficients of like-power terms, we arrive at 

A (Ur/hciy-^ 


fA(r) = - 


/b(0 = 7 


1 + wAj (Urlhc,r-m 

B {Ur/hcif^-^ 


-|l /2 


-|l /2 


(44) 


(45) 


1 + sSS) (C'r/fc, )*<'-■/">_ 

where it can be proven by construction that no such solution exists when /t 7 ^ 0 . 
Additionally, the constants in Eqs. (44)-(45) satisfy the constraint and normalization 
condition 




B 


\B\^B 

CA 


= 4 £ - 2 , 


rdr 


|A |2 (Ur/hciY^-^ + 

\B\ 

[Ur/hci)^^ ^ 


) + Am (Ur/nc,)o(‘-^i£ 



= 1 


(46) 


When the winding number £ = 0 or 1 an angular momentum term appears in only 
one of Eqs. ([^-(|^. In this case the NLDE is easier to solve and allows for algebraic forms 
even when /i 7 ^ 0. For the homogeneous case ju = 0 though the results in Eqs. (44)-(45) 
can be used to obtain the correct solution. Substituting £ = 1 into Eqs. (44)-(45) gives 
the upper and lower spinor components 

A 


fA(r) = - 


1 + ^5 ([/r/fe,)* 


11/2 


(47) 
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Figure 4. (color online) NLDE vortex radial solutions, (a) Bessel solution for i = 3; 
(b) numerical solution for 1 = 4; (c) vortex/soliton; (d) Anderson-Toulouse vortex; 
(e) ring-vortex solution for ^ = 4; (f) ring-vortex/soliton solution; (g) Mermin-Ho 
vortex; (h) half-quantum vortex. The plot in (a) is the square modulus of the complex 
asymptotic solution in Eq. (43). In each plot, the upper and lower spinor radial 
solutions are indicated in red and blue, respectively. 


fB{r) 



BjUrlhci) 
^ ((7r//iq)i 


(48) 


Equations gg-g) describe a vortex whose density peaks in the shape of a ring with 
a bright soliton (no rotation) located at its center. This solution corresponds to the 
ring-vortex/soliton from our previous work [IHj, but obtained here as a special case of 
the more general result in Eqs. (44)-(45). Note that Setting £ = 0 in Eqs. (44)-(45) 
effectively interchanges the forms for fA and fs but leads to the same solution type. 
The ring-vortex/soliton radial profiles are plotted in Fig. |^d) with the density and 
phase shown on the right of Fig. 

Addressing now the case /r 7 ^ 0 and ^ = 1 , an algebraic solution is obtained using 
an ansatz similar to that used to obtain Eqs. (44)-(45), which gives 


fA{r) 

/b(4 


y^lJ,/U{lJ,r/nci) 

[l + (Air/hQ)2]'/'’ 

[1 (yur/hQ)2]^/^ 


(49) 

(50) 


This solution describes a coreless vortex, a vortex with unit rotation in /a and a bright 
soliton in fs centered at the core of the vortex. This solution can also be obtained by 
beginning with the ansatz = tanh[ 5 '(r)], fs = sech[ 5 '(r)] which, upon substitution 
into the NLDE, may be directly integrated to give g{r) = arcsinh(;ur/hQ); applying a 
standard identity then yields Eqs. p^-(|50l). The vortex/soliton radial solution, density 
and phase are shown in Fig. |^c) and Fig. (left side), respectively. Radial plots for the 
£ = 4 ring-vortex are shown in Fig. |^c) (left side) with density and phase plots for the 
case £ = 2 shown in Fig. (left side). 
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Figure 5. (color online) Phase and density plots for the the £ = 2 ring-vortex (left) 
and the £ = 1 ring-vortex/soliton (right). (a,b) A-sublattice density and phase. (c,d) 
B-sublattice density and phase. 


6.3. Skyrmion solutions 

To obtain Skyrmion solntions we choose an ansatz of the form Ja — 3 cosip , fs = V sin(/p, 
where the parameters rj and ip are fnnctions of the radial coordinate. For backgronnd 
on skyrmions in 2-dimensions see [121 EH- Substituting these forms into Eqs. (§-(§ 
reduces the NLDE to two hrst-order nonlinear ODEs 
dp (1 -|- £) 
dr 


2 r 


-sin2(p -b 


+ cos^2(p) - ^ , 
2hci hci 


drj £ ^ 1 2 U n . ^ 

— = r7-cos2(p — Tj-cos p -\- —rj sindiy? . 
dr r r nci 


(51) 

(52) 


The first point to note is that the centripetal terms place a restriction on the behavior 


of p for r ^ 0. Equation (51) forces the condition p ^ mv/2 (n E 2.), and the only way 


to keep Eq. (52) hnite at r = 0 is to require £ = 0,1. To satisfy these conditions we find 


two possible solutions: £ = 0 for 9?(0) = 7r/2, and £ = 1 for 9?(0 ) = tt. Thus, skyrmion 
solutions exist only for one unit of angular momentum in either the upper or the lower 
two-spinor component. As we noted previously, choosing £ = 0 or 1 simply transfers a 
unit of rotation from one component to the other. We can reduce the problem of solving 


Eqs. (51)-(52) by further adding the constraint of a constant amplitude rj = C. We then 
obtain the value for C by examining the asymptotic form of the equations for r —?• cx). 
In this limit (assuming finite energy) derivatives vanish and we obtain the asymptotic 
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values 


<^(oo) = 


rUTT 


r] = ±A 


Aja/U 


4 ’ ■' y 3 + (- 1 )™ ’ 

where m G Z. Next, we combine Eqs. (0-(|5|) into one equation for Lp 


d(p 1 . h' t 

— = -sm2(p - ^ + 
dr r hci nci 


for ^ = 1, and 
dip 


dr 




fi 2/rC'„ 


hci 


+ 


he. 


1 + 


1 + 


COS^(/P 


sin^(/p 


AfJ'Cm . ^ ^ 

—-r sin4(/p 

hci J 


AflCr, 

he. 


r sin4(/p J 


(53) 


(54) 


(55) 


for £ = 0, where Cm = 1/2 for odd m, and Cm = 1/4 for even m. Equations. (51)-(52) 
allow for two types of solutions labeled by the subscript m; one solution asymptotically 
approaches 7r/4, whereas the other solution approaches 0. The Anderson-Toulouse 
solution is obtained for (/^(O) = 7r/2 and (/?(oo) = 0 [SHI EH], whereas the Mermin- 
Ho solution corresponds to the case p(oo) = n/4 in Eq. (53) [58l |60] 


Note that for 

the Anderson-Toulouse and Mermin-Ho solutions we find that r) = ^J2ii/U for the 
constant envelope factor. The radial profiles for both the Mermin-Ho and Anderson- 


Toulouse vortices are obtained by solving Eqs. (|54|)-(|55|) using a straightforward shooting 
method 


j . We have plotted the radial solutions for both types of skyrmions in Eig. |^d) 
and (g), with plots of the density and phase shown in Fig. (right side) and Fig. (left 
side). 


6.4- Half-quantum vortices 

The NLDE supports half-quantum vortices, i.e., vortices with half-integer winding 
i = 1/2. Such solutions are obtained by forming superpositions of the components 
fA and fs which solve the Mermin-Ho condition ining. The necessary requirement 
is that both spinor components approach the same nonzero value at spatial infinity far 
from the vortex core (see Fig. |^g)). The appropriate spinor combinations read 

(56) 

(57) 


/a = —ie*^/^sin(y9 — ie , 


fs = ie'^^^sinp 


ie *^/^cos 9 ?, 


where p{r) is the Mermin-Ho parameter. Note that Eqs. (56)-(57) constitute a solution 
of the full time-dependent NLDE but do not solve the time-independent case. To see 


that Eqs. (56)-(57) are associated with a fractional winding number we note that far 


from the vortex core we have p —?• 7r/4, and the wavefunction takes the form 
T = 2 [icos(d/2), sin(d/2)]. 


(58) 


From Eq. (58) we compute the geometric phase by circling the vortex core through the 


Berry phase prescription 


exp 




( 69 ) 
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Figure 6. (color online) Phase and density plots for the ^ = 1 vortex/soliton (left) 
and i = 1 Anderson-Toulouse skyrmion (right). (a,b) A sublattice density and phase. 
(c,d) B sublattice density and phase. 



Figure 7. (color online) Phase and density plots for the 1 = 1 Mermin-Ho skyrmion 
(left) and half-quantum vortex solutions (right). (a,b) A sublattice density and phase. 
(c,d) B sublattice density and phase. 


where cfs is the Berry phase and 9 is the polar angle. The vortex wavefnnction 
transforms as a Dirac spinor under spatial rotations accumulating a phase factor 
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exp (—zcr20/2), where cr^ is the third Pauli matrix. Working out the exponential in 


Eq. (59) gives 

/*27r 




'0 


_ gi6>/2 cos(51/2) , e sin(6*/2) 


d 


ei0/2 sin(d/2) 


^ ^-ie /2 I i e cos(d/2) 


r*27r 


= / - - 


= —ZTT . 


(60) 


Radial profiles for half-quantum vortex spinor components are plotted in Fig.j^h) where 
we have taken the polar angle equal to zero. Density and phase plots are shown in Fig. 
(right) with general properties of vortices tabulated in Table 


Vortex type 

Winding 


Analytic form of T(rJ 

) 

Asymptotic 

behavior 

Vortex/soliton 

i = l 


1 ^iO {r/ro) 

T 

\iPa{oo)\ = 1, 


y^l+(r/ro)2’ \/l+(r/ro)2 


|k’B(oo)|=0 

Ring-vortex/soliton 

i = l 


^ 1 ^ie (r/ro) 

T 

|■0.4(oo)| =0, 


^l+(r/ro)^’ yfl+{r/roY 


|V’b(oo)|=0 

Anderson-Toulouse skyrmion 

i = l 

[ 

z cos(/?(r/ro), e^^sin(/9(r/ro)]^ 

(y9(oo) =0 

Mermin-Ho skyrmion 

i = l 

cos(/?(r/ro), e^^sin(/?(r/ro)] 

(^(cxd) =7r/4 

Half-quantum vortex 

i = l 


[zcos6>/2, sin6>/2]^ 


|«'(oo)| = l 

Ring-vortex 

^ = 2,3,4,... 

Aj(e-i)0_ 

(r/ro)^ ^ Me {r/roY^ ^ 

^ IV'a(oo)|=0, 

'y/l+(r/ro)®^^ 1/2)’ 1/2) 

|k’B(oo)|=0 

General topological vortex 

£ = 2,3,4,... 

Numerical shooting method 

|■0A(oo)| = 1, 






|V’b(oo)| = 1 


Table 2. Vortex solutions of the NLDE. Solutions are described by their phase 
winding, closed-form expression, and asymptotic properties. The Mermin-Ho, 
half-quantum vortex, and general topological solutions have conserved topological 
charge associated with azimuthal phase winding classified according to elements of 
the fundamental group 7ri(5'^) = Z. In all other cases, the spinor components 
asymptotically approach either a saddle point or a local maximum of the effective 
potential. Note that ro is the length scale associated with the chemical potential or 
the interaction strength depending on the particular solution. 


7. Vortices in the semiclassical reduction to nonlinear Schrodinger equation 
with correction terms 

We now study the connection of NLDE vortices to their nonlinear Schrodinger 
counterpart. Working again from the ansatz '0^(^)(r,t) = where we 

have factored out the dominant energy contribution from the spinor functions, and 
substituting into the NLDE gives 

Iiva + to/1,1 = -iVvB + U\vb\^Vb , 

^ivb + iVB,t = -iVvA + U\vb\'^vb , 


(61) 

(62) 
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where we have abbreviated the differential operators by using V = d^-l- idy. For clarity 
of notation, we have set h = q = 1 and used the abbreviated subscript notation to 
indicate differentiation. We then make the following approximations 

1. \vb\^ ^ 

2. \vB,t\ -c \va,x\ 

3. \U\ -C p. 

With these approximations, Eq. 

i 


can be solved for vb to give 


vb 




-{Vva) 


1 + —\va\‘ 

fl 


Substituting Eq. (6^ back into Eq. 
i 


ivA,t = 


/i 


'-{Vva) + 


U\va\'^va - iiva , 


(63) 

(64) 

(65) 

( 66 ) 

(67) 

( 68 ) 

(69) 


= - i (| P |^ J > j ) (l + —- Vvva){V‘\va\^) + U\va\^va - fiVA , 

= -i(V^wj) (l + —|WA|d - tt{'DVA){'D‘\VA\'^) + UIvaI'^va - flVA ■ 

II V 

Note that equality from this point on must be understood within approximations 
Eqs. (63)-(65). Einally we obtain 

iVAf = -tv^VA + UlvAf^VA - i^va “ ^ + (PojXP*IwaH] ■ 


(70) 


Equation (70) is the nonlinear Schrodinger equation with derivative correction terms 
and attractive constant background potential —/r. Note that the effective particle mass 
is /i/2. The correction terms in Eq. (70) are smaller than the first three terms on the 
right hand side due to the extra factor oi U/(i and powers of derivatives of va- This 
tells us that localized solutions will differ from standard NLS form only in regions where 
gradients are steep, i.e., near the core of a vortex. 

Next, we look for vortex solutions of Eq. (70) by choosing the stationary, 
axisymmetric factorized form 


VA{r,e,t) = Ce-^‘^*V^%{r), 
where (7 G R and n G Z are constants, and such that 


lim, 


.\va\=C 


lim, 


|j>(r)| = 1 . lim^oIxAl = 0. 


(71) 

(72) 


In plane-polar coordinates, Eq. (70) becomes 

1/^2 I Id 

ji 86“^ r dr 


^VA,t 


Va + {U\va\‘^ - fr)vA 


U 


{ d^ I d^ 




1 d 


u 

~2 




JO 


d .Id 
r d6 


Va 


■,-iO 


d .1 a , 2 


(73) 
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Inserting Eq. 0 gives 

1 


OJV = —- 




dr^ 
f 


--^CV + 


n 


n 




2^ ^ V dr‘^ 


r dr 
r dr 


V + {UC^v^ — iJi)v 


p? 


d n 
dr r 


d n\ / d n 
dr r / \dr r 


(74) 


Canceling some terms and using condensed notation, we have 

1 / .2 


UJV = - 

p 

A 


n“ 1 


n 


1 


Vrr - 7;V + -Vr 

fytZi 


More simplifying leads to 


1 r 1 

— s Uj-j. H- Vr + 

p [ r 

= ucA^ - 

p^ 


p{p-\- (jj) - 


+ {UCA^ - p)v 


U n2 ( ^ 

- -C V (Vr - V 

jji^ \ r 


n 


2Vr 


(75) 


A 1 

^rr 


r 


r 


o ^ ^2 (2 ^ 

- 2—C V (v^ - vvj 

jji^ \ r 


(76) 


This can be further simplified by making the coordinate change ^ p{p-\- lv) r 


(p + U)) 


Ai + + 


A 

e 


u 


= UC^v'^ - {p + u) “ To 


p 


n 


V + - 


e e 


A 


—2— {p + Lu) C^v • 

Finally, dividing through by (p + lj) gives 


(77) 


1 

Ai + + 


P 


n 

e 

n 

A - 


UC^ 3 ^^2 2 

7 -r - 

yp + UJ) p 


n 


1 


V(( - -^v + -Vi 


(78) 


= 0 . 


(79) 


There are several limits of this equation that are interesting: 

(i) For small r we expect that n —?• 0. In this limit the nonlinear terms on the right 
hand side become negligible so that 

1 / A 

This is Bessel’s equation and the solutions are the well known Bessel functions JnA) 
and Yn{x). The ones we are interested in are the Bessel functions of the first kind, 
Jn{x), since these are regular at the origin. Of these, we further restrict solutions 
to those for which n > 0 since these are zero at the origin. 

(ii) We seek solutions with constant square modulus for large r. In this limit, all 
derivatives may be set to zero in Eq. (|7^. This allows us to solve for C in terms 


of the constants U, p, and uj. For large ^ (large r), n —?• 1, and Eq. (78) reduces to 
1 = UC^Kp + C(j), or simply C = ^{p-\-uj) /U. 
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(iii) For U ^ 0, we should retrieve the standard Schrodinger result for the radial part of 
the wavefunction. This is indeed the case since, in this limit, all terms on the right 
hand side of Eq. (78) vanish and we are left with Bessel’s equation as expected. 


8. Vortices by the method of numerical shooting 

In general, vortex profiles for arbitrary phase winding i. and chemical potential fi may be 
obtained using a numerical shooting method as described in Refs. |13] and [HH!- These 
include topological solutions whose tails do not vanish asymptotically, as well as non- 
topological vortices where one or both spinor functions have asymptotically vanishing 
profiles. To proceed we first convert the NLDE to dimensionless form by introducing 
the rescaled radial coordinate and spinor components 

X = lar/(hci ), riA = y/uJJifA , Vb = , (80) 

so that Eqs. (|8|)-(@) become 

~ = -Va{x) , ( 81 ) 

- ~ Mx)fvBix) = -VbU) ■ (82) 

Here the dependence of the solution on the choice of angular quantum number i is 
implied. The asymptotic form of these equations shows that moduli for convergent 
solutions \riA(B)\ can approach either 0 or 1 for large %• To study the analytic structure 
as well as a practical starting point for application of numerical methods, we expand 
the solution in a power series in % 


oo oo 

Va (x) = = ■ (83) 

j=0 j=0 


Substituting these forms into the NLDE gives relations for the 

expansion coefficients 

O 

1 

II 

CO o 

1 

+ 
T—1 

(84) 

(2 “h — SdiQjQ — — Ci\ , 

(85) 

(3 + £)bs — 3a2al — = —a2 , 

(86) 

(2 — £)ai + 6q = 6o , 

(87) 

(3 - i)a2 + Sbibl = bi , 

(88) 

(4 — £)cis + 362^0 360&1 = b2 , 

(89) 


Choosing a particular value for i determines the values of both Oq and bo and we find 
that there is only one independent parameter. Thus, for a given value of £ a vortex 
is found by tuning a^_i towards a critical value As examples, we have found 
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Ur/hci Ur/hci Ur/hci 

Figure 8. (color online) Radial profiles for i > 1 vortices by numerical shooting, (a) 
Profile for ^ = 2. (b) Profile for ^ = 3. (c) Profile for = 4. 


vortices for the three lowest £ values for which both spinor components have 
rotation 

nonzero 

^vortex ^ Q ^ 

£ = 2, 

(90) 

^vortex ^ 

£ = 3, 

(91) 

ojortex ^ 0.0240267... , 

£ = 4. 

(92) 


The associated radial profiles for these solutions are plotted in Fig. In Fig. we 
illustrate convergence for both spinor components for the £ = 2 vortex. Convergence to 
the vortex profile is seen by overlaying an excitation of the vortex (a) and the free-particle 
Bessel-like solution (b). Note that there are two types of radially excited solutions 
characterized by the property that |?]u,(b)| oscillate around 1 or around 0. 



Figure 9. (color online) Convergence of numerical vortex for £ = 2. (a) For 

oi > the solution overshoots to a radial excited state of the vortex. Note that 

such oscillating solutions converge to unity far from the core, (b) For Ui < the 

solution undershoots and converges to the linear solution Bessel functions. The solid 
blue and dashed red plots are the A and B sublattice radial wavefunctions, respectively. 
The solid black and dashed black plots are the exact solutions for the A and B sublattice 
radial wavefunctions, respectively. 


We may elaborate further on the oscillating regimes on either side of the fiat vortex 
profile by examining the solution for £ = 1. The £ = 1 state corresponds to one unit of 
rotation in ipB and no rotation in if a, in which case a trivial constant solution exists, i.e., 
= 1 and r^B = 0, for0>x>oo. We address this solution as it offers a clear look at 
higher radial excitations. To see this we solve Eqs. (81)-(82) by numerical shooting for 
£ = 1. The resulting spin component profiles are plotted in Fig.[T^a)-(c) where we have 
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also included plots of the associated total densities + in Figs.[lO)[^d)-(f). 

Radial excited states appear as we tune the initial value for riA, or equivalently Oq in the 
Taylor expansion, away from unity. The spatially constant solution obtained for oq = 1 
can be thought of as a boundary between two oscillating regimes: one characterized 
by strong nonlinearity with oscillations around r/A = 1, riB = 1 in Fig. Htb) ; the 
other associated with weak nonlinearity with oscillations around tja = f]B = ^ shown in 
Fig. [T^c). When tuning Oq upward from unity, the total density in Fig. [I0|^e) shows 
the onset of inward movement of bright vortex rings over a nonzero background. For 
this plot we have chosen Uq = 1 + 10“^. Here we see that the solution overshoots to an 
excited state of the vortex and enters a regime where the nonlinearity dominates the 
kinetic energy. In contrast, for the choice oq = 1 — 10“^ Figs. [To|c) and (f) show the 
onset of radial oscillations which result from a dominant kinetic energy in the equations 
of motion. In this case the total density exhibits decaying oscillations towards a zero 
background density. Here, undershooting the flat solution results in solutions which 
approach the Bessel function form as seen in Fig. [Io|c). This is what we expect to 
see for weak nonlinearity. This illustration highlights the property of the Dirac kinetic 
terms which act by confining the solution in the presence of strong nonlinearity 



Ur/hci Ur/hci Ur/hci 

Figure 10. (color online) Radial excitations for i = 1. (a) Globally flat solution, (b) 
Strong nonlinearity, (c) Weak nonlinearity. 

So far in this section we have shown how tuning the initial condition above some 
critical value forces the solution into the strongly nonlinear regime, whereas below this 
value leads to solutions in the weakly nonlinear regime. In our analysis we held the 
chemical potential and interaction fixed, fi = U = 1. We now study the effect of tuning 
the chemical potential for p > 1 —)■ /r = 0, while maintaining the initial value in the 
shooting process and the interaction fixed. For arbitrary phase winding we will see 
a progression from Bessel solutions through a topological vortex, then finally ending 
in a single bright ring-vortex with vanishing tail. This solution resembles those found 
in attractive BECs and also in optics. Although we demonstrate such solutions by 
analytical methods in the next section, the numerical approach allows ns to demonstrate 
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the transition from the free-particle limit, /a ^ U (weak nonlinearity), to the vortex 
limit, fj, U (strong nonlinearity). To show this transition we choose a particular 
rescaling of the NLDE such that the limit /r ^ 0 makes sense, i.e., 

X = U‘^r/{hciia ), r]A = y/pjuJ a , Vb = • (93) 

The resulting dimensionless NLDE is 


- a 


a. 


£ 

'y H - 

^ X 

1 -1 

h- 

X 


Vb{x) + \’lA(x)f VAix) = fiilAix), 
Va(x) + l>7B(X')t>7B(x) = naix) ■ 


(94) 

(95) 


Starting with the case £ = 1, we fix Uq = 1 and 6o = 0 in the Taylor expansions, 
Eq. (83), and tune ft = pi/U toward zero starting from jl > 1. The progression 
of this solution as /i ^ 0 is depicted in the sequence of plots in Fig. [T^ As Ji is 


reduced towards 1, Fig. [TT|(a)-(c), Bessel-like oscillations about riA{B) = 0 are pushed 
out towards large X; completely flattening out the solution at p = 1 in Fig. [II|(d). As 
we continue decreasing jl towards 0, oscillations about riA(B) = 1 move inward from 
large x and finally flatten out leaving only a coreless ring-vortex centered at x = 0 
shown in Fig. [IT|(i). The precise values of the renormalized chemical potential are 


jl = 1.22, 1.09, 1.0005, 1, 0.9995, 0.775, 0.55, 0.316, 0.1, for the panels (a)-(i) inFig.jn 
respectively. 

Analogous progressions are seen for arbitrary £ values culminating in a single ring- 
vortex centered at x = 0 when jl = 0. The qualitative difference between the solution 
for £ = 1 and those for £ > 1 is that in the later case both components must vanish at the 
origin. In terms of numerics this requires taking bo = oq — 0, while specifying the first 
derivative of tja at the origin. The progression is plotted for the case £ = 2 in Fig. [T^ 
Note that for each value of £ there are two distinct types of asymptotically fiat vortices, 
as displayed in panels (d) and (i) of Fig. 12, The specific values of the renormalized 
chemical potential are jl = 1.3, 1.21, 1.2047, 1.204267325, 1.203, 0.85, 0.55, 0.32, 0.1, 


for panels (a)-(i) of Fig. 12 


The NLDE allows for vortex solutions which satisfy an additional symmetry given 
by the constraint ItjaI"^ + \vb\‘^ = 1- Solutions which satisfy this constraint are the 
vortex/soliton or coreless vortex, Anderson-Toulouse and Mermin-Ho skyrmions. Such 
solutions may be obtained numerically by shooting backwards from large x towards 
X = 0. In particular, for the Mermin-Ho solution we integrate backwards starting with 
the asymptotic boundary condition tja = cos(7r/4) -|- 10“^, riB = sin(7r/4) — 10“^. Here 
k is a. parameter tuned to give the desired values of functions at the origin, analogous 
to ao for the forward shooting. The Anderson-Toulouse vortex is obtained similarly but 
with the boundary condition rjA = cos(7r/2) -|- 10“^, tjb = sin(7r/2) — 10“^. The half- 
quantum vortex is obtained by forming linear combinations of the numerical Mermin-Ho 
components, as previously discussed. 
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Figure 11 . (color online) Progression of radial profiles for vortiees with phase winding 
£ = 1. (a)-(c) Solutions where the renormalized chemical potential satisfies the 

condition /i > 1. (d) Solution for /i = 1. (e)-(i) Solutions for /i < 1. The sequence of 
plots depicts the transition from the weakly-interacting Bessel-like solution in (a) for 
which /i > t/, to the strongly nonlinear case for which p < U in {i). Note that ft = 1 is 
the boundary between solutions which oscillate around 0 and solutions which oscillate 
around 1. The results in these plots show that general ring-vortex solutions may be 
obtained by starting from excited states such as in (e), then reducing p towards zero. 


9. Discrete Spectra in a Harmonic Trap 

We now extend our numerical studies from Sec. [8] to the treatment of vortices bound 
within a weak harmonic potential similar to the analysis in our previous work on solitons 
(a brief summery of these results can be found in [H], presented here in full detail). 
Physically this accounts for an additional external harmonic trapping potential present 
in real experiments. We follow a similar procedure as before in preparing the NLDE 
for numerical analysis by first converting to a dimensionless form. Here we examine 
the same solutions as in Sec. but with tight confinement assumed in one of three 
spatial dimensions and with additional weak harmonic confinement in the remaining 
two directions. The oscillator frequencies thus satisfy oj = ujx^ ujy. Equations. 
i§ are already defined for a quasi-2D system, as the z-dependence has been integrated 
out and parameters are normalized accordingly. Thus, we require the harmonic potential 
to be dependent only on the planar directions x and y. We then introduce the planar- 
symmetric harmonic potential V{r) = {1/2)MixP‘{x^ + = (l/2)Ma;^r^. Next, we 

choose a dimensional rescaling of the NLDE appropriate to the harmonic oscillator. 
We divide through by the harmonic oscillator energy huj and define the dimensionless 



















Figure 12. (color online) Progression of radial profiles for vortiees with phase winding 
i = 2. (a) /i = 1.3. (b) ft = 1.21. (c) ft = 1.2047. (d) ft = 1.204267325. (e) jl = 1.203. 
(f) jl = 0.85. (g) jl = 0.55. (h) jl = 0.32. (i) jl = 0.1. The topological vortex solution 
is shown in (d) and isolated ring-vortex in (i). 


variable and spinor components in terms of this energy scale 
X = hujr/{hci) , r}A(B) = y/U/hhJ fA(B) ■ 

This transforms the NLDE to 

- + ^Vb{x) + \VA{x)fvA{x) + Qx^rjAix) ^ I^VAix), 




l-£ 

X 


Va{x) + \VB{x)fvBix) + Qx^Vb{x) = fiVBix) 


where the two dimensionless parameters in our equations are 

Q 


Mcf 


2 fkj 




JL 

hu) 


(96) 

(97) 

(98) 

(99) 


Solving Eqs. (97)-(98) subject to the radial boundary conditions consistent with the trap 
potential leads to radially quantized solutions labelled by integer quantum numbers. 
Fig. [T^ shows the first six quantized solutions for the i = 2 topological vortex in panels 
(a)-(f), respectively. 

To connect our results to experiment we compute the discrete eigenvalue spectra 
for our radially quantized solutions which relate the chemical potential for a particular 
solution to the strength of the nonlinearity. The spectrum is computed by first defining 
the normalization condition as an integral over the dimensionless spinor components 


/ 


xix{\<iA{x)? + \'nB{x)?) = JV, 


( 100 ) 






















Vortex solutions and spectra in a weak harmonic trap 


30 


fA 

fB^ 

0 

-0.5 

fA 

1 

/b 

0 

-1 

02468 10 02468 10 0 246 8 10 

Ur/hci Ur/hci Ur/hci 

Figure 13. (color online) Numerical shooting for quantized vortices in a quasi- 
2D harmonic potential, (a)-(f) First six quantized states for ^ = 2. The extreme 
vortex limit (a) is characterized by a dominant nonlinearity with the spinor functions 
flattening out over a greater portion of the domain with a weaker contribution from 
the derivate terms. In contrast, in the free-particle limit (f) spinor functions resemble 
Bessel oscillations. 







where V on the right hand side encapsulates key lattice, condensate, and trap 
information through the relation 

VShwNU 


w = 




( 101 ) 


The calculation proceeds by fixing the value of Q in Eqs. (97)-(98) and varying jl, while 
computing the norm W for each value of fl by integrating Eq. (100). This procedure gives 
the paired values Eixing the total particle number N gives a relation between 

the chemical potential fi and the interaction U. The values for the free parameter in the 
Taylor expansion, uq, normalization A/", and corresponding chemical potential fl = [i/hw, 
are tabulated in Table for the lowest radial i = 2 mode shown in Eig. [T^a). Plots of 
this data along with the spectra for other solutions are shown in Fig. [T^ We have taken 
Q = 0.001 for all of our calculations. The spectrum for a fixed quantized mode tracks 
the flow of the chemical potential as one tunes the system between the free-particle and 
strongly nonlinear limits. In the former limit radial profiles for the lowest excitations 
resemble decaying Bessel functions as in harmonically trapped massless Dirac particles, 
whereas in the latter limit the nonlinearity dominates and the lowest trapped solutions 
flatten, consistent with the Thomas-Fermi regime. For higher quantized modes radial 
derivatives force the solution to maintain the Bessel-like form even as one tunes the 
parameters into the strongly nonlinear regime. Convergence of our solutions is provided 
in Appendix A. 
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Free parameter uq 

Normalization Af 

Chemical potential jl 

0.00000003 

1.9679 X 10“® 

1.8 

0.00000365 

0.29560 

2 

0.0000 1018 

1.85630 

3 

0.00001532 51 

3.63288 

4 

0.0000208947 

5.87547 

5 

0.00002708 8395 

8.62326 

6 

0.000033870293 

11.72273 

7 

0.00004118 586175 

15.54209 

8 

0.000048991422 392 

20.12707 

9 

0.0000 5725 27244133 

24.87695 

10 

0.000065942040680759 

30.25986 

11 

0.000075036 2691057811 

35.51507 

12 

0.0000 84515722 6783 7827 3795 

37.30714 

13 


Table 3. Numbers for computing spectra of i = 2 topological vortex ground state 
solution. For a fixed value of the chemical potential fr fhe free parameter oq is tuned 
until the desired radially quantized state is reached, for which one then computes the 
normalization Af by Eq. (100). As fi increases into the Thomas-Fermi regime the 
dependence of the solution on oq becomes more sensitive, requiring a greater degree of 
tuned accuracy as shown in the column on the left. 





Figure 14. (color online) Spectra for relativistic vortices confined in a harmonic 
potential.(a) Vortex/soliton (black curve), Anderson-Toulouse skyrmion (red), Mermin- 
Ho skyrmion (blue), and half-quantum vortex (green), (b) General topological vortices 
for = 2,3,4 (black, red, blue), (c) Radial ground state and first two excited states 
of the vortex without skyrmion symmetry (black, red, blue). In each figure, the 
renormalized chemical potential is plotted as a function of the normalization. 


10. Conclusion 

In this article we have solved the NLDE by a variety of different methods for both the 
idealized case of zero trapping potential and in the harmonic case. The latter context 
provides the opportunity to study spatial quantization of vortex solutions in the radial 
dimension. A preliminary asymptotic solution using Bessel function expansions provides 
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insight into the strncture of the NLDE itself. Algebraic solutions were then obtained by 
considering the case of zero chemical potential. For these solutions the derivative and 
nonlinear terms are perfectly balanced leading to bright vortex rings over a zero-density 
background, a direct consequence of the Dirac operator. More generally, for nonzero 
chemical potential, we found that analytical solutions are only possible when one unit of 
winding is considered. In this case we used a numerical approach to obtain vortices with 
arbitrarily large winding number. A combination of numerical and analytical techniques 
yields skyrmion and half-quantum vortices, i.e., textures. 

Having obtained our solutions, we computed their discrete spectra in the presence of 
a weak harmonic potential. This gives us the low-temperature /r versus U landscape for 
relativistic vortices, where p and U are the chemical potential and lattice renormalized 
particle interaction, respectively. For example, for finite pL we found that a series of phase 
transitions occur as U is tuned from zero upward: we encounter a Mermin-Ho skyrmion 
transitioning into a half-quantum vortex, followed by the Anderson-Toulouse skyrmion, 
then finally into a vortex/soliton, i.e., a bright soliton at the core of a singly-wound 
vortex. 

Some of our solutions are similar to those obtained in spinor BFCs, skyrmions in 
particular. However, presence of the Dirac operator contrasts heavily with the Laplacian 
case. This difference is most obvious in our bright ring-vortex solutions. Components of 
these solutions resemble the bright vortex which occurs in attractive BFCs, but in our 
case the confining (focusing) regime of the Dirac operator is responsible for the effect in 
spite of repulsive atomic interactions. Furthermore, we should emphasize the significant 
distinction between our results and similar confinement effects in spinor BFCs and in 
traditional models such as Thirring and Gross-Neveu. In each of these other theories 
the presence of confinement relies on attractive interactions in addition to a mass term, 
whereas we consider strictly the repulsive case with zero mass. 

The interdisciplinary nature of our work suggests several future research directions. 
For instance, the form of our solutions implies possible mappings to finite energy 
solutions to classical gauge field equations. In particular, there is a deep connection 
between the NLDE and Chern-Simons terms relevant to quantum Hall fluids and more 
generally to relativistic field theories [HU ESI ES] • To see this connection consider that 
Eq. 0 resembles the low-energy effective theory for massless Dirac fermions interacting 
through a gauge field, where the latter has been absorbed into the local fermion 
contact terms. The Dirac term by itself is symmetric under global phase and chiral 
rotations which are made local through the addition of a gauge field. Quantization 
results in the well known axial anomaly which one finds to have exactly the Chern- 
Simons structure |67l EH ES] . Thus, quantized theories of interacting fermions generally 
result in couplings to Chern-Simons terms. In our case we have focused only on the 
fermion part of the argument which naturally retains imprints of the omitted Chern- 
Simons terms. Another argument for the NLDE/Chern-Simons connection hinges 
on the well established duality between the Thirring model in (2-|-l)-dimensions and 
Maxwell-Chern-Simons theory. The mapping is arrived at through bosonization, i.e.. 
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reformulation in terms of paired particle and antiparticle field operators effective at 
strong coupling or low energy. These points merit further inquiry into potentially 
fruitful questions. Finally, in light of the successful analogs to date connecting condensed 
matter and nonlinear optics we expect that our results should be reproducible within 
an optics setting. Another obvious direction is to go deeper into the mathematics of 
nonlinear partial differential equations. The vast amount of work in this area provides 
an array of solution and classification methods which may be used to fully understand 
critical bounds for well-posedness as well as general solutions to the full NLDE with 
time dependence. 
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Appendix A. Convergence of numerical vortex solutions 


We demonstrate convergence of solutions in the harmonic trap for three of the 1 = 2 
topological vortices associated with the black curve in Fig. [T^b). Radial profiles for 
the chosen solutions are shown in Figs. |Al[ a)-(c). The corresponding values of the 
chemical potential are /i = 4, /r = 7, and /r = 10. These values interpolate between 
the free-particle and strongly nonlinear limits (small to large /i values). These solutions 
were obtained by finite differencing using a shooting method to tune the precision of the 
initial value of ipA such that <C 1 to pick out the ground state. For convergence at 
a single radial point, we compute the value of the solution at the dimensionless radius 
Xi = Ti/^Dij-sc = 10 for several values of the grid size N = 10^, 10^, 10^, 10^, 10®. We use 
the error formula which depends on the dimensionless radius and number of grid points 


N) = 


+^(xOn(B)) 


(A.l) 


where in the symbol 'ip{Xi)\B) subscript A{B) denotes the sublattice excitation, 
Xi denotes the element in the discretized dimensionless radial coordinate, and the 


superscript N denotes the number of grid points used in the calculation. In Figs. Al(d) 
(f), we have plotted log^g N)| versus log^^p 

Figs.|ATKa)-(c). 


N, for the solutions shown in 
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